home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / Mathematics / Notebooks / COSY_PAK / SignalProcessing / Analog / InvLaPlace.m next >
Encoding:
Text File  |  1992-06-08  |  39.5 KB  |  1,188 lines
  1. (*  :Title:    Inverse Laplace Transform  *)
  2.  
  3. (*  :Authors:    Wally McClure, Brian Evans, James McClellan  *)
  4.  
  5. (*
  6.     :Summary:    Multidimensional, bilateral inverse Laplace transform
  7.         for signal processing expressions.
  8.  *)
  9.  
  10. (*  :Context:    SignalProcessing`Analog`InvLaPlace`  *)
  11.  
  12. (*  :PackageVersion:  2.6    *)
  13.  
  14. (*
  15.     :Copyright:    Copyright 1990-1991 by Brian L. Evans
  16.         Georgia Tech Research Corporation
  17.  
  18.     Permission to use, copy, modify, and distribute this software
  19.     and its documentation for any purpose and without fee is
  20.     hereby granted, provided that the above copyright notice
  21.     appear in all copies and that both that copyright notice and
  22.     this permission notice appear in supporting documentation,
  23.     and that the name of the Georgia Tech Research Corporation,
  24.     Georgia Tech, or Georgia Institute of Technology not be used
  25.     in advertising or publicity pertaining to distribution of the
  26.     software without specific, written prior permission.  Georgia
  27.     Tech makes no representations about the suitability of this
  28.     software for any purpose.  It is provided "as is" without
  29.     express or implied warranty.
  30.  *)
  31.  
  32. (*
  33.     :History:    date started -- April 25, 1990 (adapted from "InvZTransform.m")
  34.         added Dialogue ability -- September, 1990
  35.         handles real-valued polynomials -- March, 1991
  36.  *)
  37.  
  38. (*  :Keywords:    Laplace transform, region of convergence    *)
  39.  
  40. (*  :Source:    Operational Mathematics by Churchill  *)
  41.  
  42. (*
  43.     :Warning:    Works only on one dimensional and separable
  44.         multidimensional transforms.
  45. *)
  46.  
  47. (*  :Mathematica Version:  1.2 or 2.0  *)
  48.  
  49. (*  :Limitation:  *)
  50.  
  51. (*
  52.     :Discussion:  1-D rules base has 70 rules:
  53.  
  54.             I.   rational transforms       26 rules *
  55.             II.  non-rational transforms   28 rules
  56.             III. transform properties      10 rules *
  57.             IV.  transform SP structures    0 rules
  58.             V.   transform strategies       6 rules
  59.  
  60.           * a partial fractions strategy rule exists in section
  61.  
  62.     Unlike some symbolic transform packages, these rules are maintained
  63.     in  a  list  called  InvLaPlaceRules.   This  was necessary because
  64.     Mathematica reordered  these  rules  if they were coded as a set of
  65.     transformation functions.  Another benefit of the list form is that
  66.     it allowed more control over how rules are applied.
  67.  
  68.     The  InvLaPlace  object (function) first  calls the interface rules
  69.     (InvLaPlaceInterfaceRules), which  handle  default arguments, error
  70.     messages, and multidimensional transforms.  These rules simply call
  71.     MyInvLaPlace  once  per  dimension  of the transform.  This inverse
  72.     transform is biased toward causal systems.   However, by specifying
  73.     the proper  region of convergence (ROC),  anti-causal functions can
  74.     be obtained (see InvLaPlace::usage).
  75.  
  76.     The  driver  for  the one-dimensional rule base is the six argument
  77.     version  of  myinvlaplace.    myinvlaplace[X, s, t, rm, rp, st, op]
  78.     either returns  the  inverse transform as a function of time (t) or
  79.     an approximation to the actual inverse transform.  If the rule base
  80.     can do neither,  which  should  not  happen, the interface function
  81.     cleanup will report any errors.  Arguments of myinvlaplace:
  82.  
  83.     X  laplace-transform function        rm  Rminus of ROC
  84.     s  laplace-transform variable        rp  Rplus of ROC
  85.     t  "time"-domain variable        st  local state semaphores
  86.  
  87.     At  each  step  in the transform rule base,  the current expression
  88.     has a local state  associated  with  it.   This state consists of a
  89.     list of six boolean values.   Each boolean value is associated with
  90.     a strategy.  If an element is True, then that strategy has not been
  91.     tried yet; if False, then that strategy has already been tried, and
  92.     it will not be tried again.   Thus,  local state is used to prevent
  93.     infinite loops which would be caused by the  repetitive application
  94.     of certain strategy rules.  See the section below called  S T A T E
  95.     D E F I N I T I O N  and see section VI of the rules.
  96.  *)
  97.  
  98. (*  :Functions:     InvLaPlace  *)
  99.  
  100.  
  101. If [ TrueQ[ $VersionNumber >= 2.0 ],
  102.      Off[ General::spell ];
  103.      Off[ General::spell1 ] ];
  104.  
  105.  
  106. (*  B E G I N     P A C K A G E  *)
  107.  
  108. BeginPackage [    "SignalProcessing`Analog`InvLaPlace`",
  109.         "SignalProcessing`Analog`LSupport`",
  110.         "SignalProcessing`Support`TransSupport`",
  111.         "SignalProcessing`Support`SigProc`",
  112.         "SignalProcessing`Support`ROC`",
  113.         "SignalProcessing`Support`SupCode`",
  114.         "SignalProcessing`Support`FilterSupport`" ]
  115.  
  116. (*  B E G I N     P A C K A G E  *)
  117.         
  118. (*  U S A G E     I N F O R M A T I O N  *)
  119.  
  120. InvLaPlace::usage =
  121.     "InvLaPlace[f, s] and InvLaPlace[f, s, t] gives the multidimensional \
  122.     bilateral inverse Laplace transform of f. \
  123.     A region of convergence can be specified by using \
  124.     InvLaPlace[{f, rm, rp}, s, t], where rm is R- and rp is R+ \
  125.     in the region (strip) of convergence: R- < Re(s) < R+. \
  126.     Note that InvLaPlaceTransform is an alias for InvLaPlace."
  127.  
  128. (*  E N D     U S A G E     I N F O R M A T I O N  *)
  129.  
  130.  
  131. Begin[ "`Private`" ]
  132.  
  133.  
  134. (*  U S E R     I N T E R F A C E  *)
  135.  
  136. (*  InvL operator  *)
  137. Unprotect[InvL]
  138. InvL/: TheFunction[ InvL[ s_, t_ ][ f_ ] ] := InvLaPlace[ f, s, t ]
  139. Protect[InvL]
  140.  
  141. (*  InvLaPlace  *)
  142. Options[InvLaPlace] :=
  143.     { Apart -> Rational, Definition -> False,
  144.       Dialogue -> True, Simplify -> True,
  145.       Terms -> 10, TransformLookup -> {} }
  146.  
  147. (*  bridge from outside world to interface to rule base  *)
  148. InvLaPlace[ x_ ] :=
  149.     invlaplacedriver[ Options[InvLaPlace], x ]
  150. InvLaPlace[ x_, s_ ] :=
  151.     invlaplacedriver[ Options[InvLaPlace], x, s ]
  152. InvLaPlace[ x_, s_, t_ ] :=
  153.     invlaplacedriver[ Options[InvLaPlace], x, s, t ]
  154. InvLaPlace[ x_, s_, t_, op__ ] :=
  155.     invlaplacedriver[ ToList[op] ~Join~ Options[InvLaPlace], x, s, t ]
  156.  
  157. (*  driver for interface to rule base  *)
  158. invlaplacedriver[ op_, rest__ ] :=
  159.     Block [    {},
  160.         tlist = {};            (* global variable *)
  161.         cleanup[ invlaplace[op, rest],
  162.              Replace[Dialogue, op],
  163.              TrueQ[Replace[Simplify, op]] ] ]
  164.  
  165. (*  interface to rule base   *)
  166. invlaplace[ op_, args__ ] :=
  167.     Block [ {},
  168.         Replace[ linverse[op, args], InvLaPlaceInterfaceRules ] ]
  169.  
  170. InvLaPlaceInterfaceRules = {
  171.     linverse[ op_, e_ ] :>
  172.         invlaplace[ op, e, LVariables[e] ] /; LTransformQ[e],
  173.     linverse[ op_, e_ ] :>
  174.         Message[ InvLaPlace::novariables ],
  175.  
  176.     linverse[ op_, e_, s_ ] :>
  177.         invlaplace[ op, e, s, DummyVariables[Length[s], Global`t] ] /;
  178.         validvarQ[s],
  179.  
  180.     linverse[ op_, e_, s_, rest___ ] :>
  181.         Message[ Transform::badvar, "s", InvLaPlace, s ] /;
  182.         ! validvarQ[s],
  183.     linverse[ op_, e_, s_, t_, rest___ ] :>
  184.         Message[ Transform::badvar, "time", InvLaPlace, t ] /;
  185.         ! validvarQ[t],
  186.  
  187.     linverse[ op_, e_List, s_, t_ ] :>
  188.         invlaplace[ op, MakeLObject[e, s], s, t ],
  189.     linverse[ op_, e_, s_Symbol, t_ ] :>
  190.         MyInvLaPlace[ e, s, t, op ],
  191.     linverse[ op_, e_, s_List, t_ ] :>
  192.         multiDInvLaPlace[ e, s, t, op ]
  193. }
  194.  
  195. cleanup[ trans_, dialogue_, simplify_ ] :=
  196.     Block [    {adjtrans, lt, sdialogue},
  197.         sdialogue = SameQ[dialogue, All];
  198.         lt = If [ simplify,
  199.               SPSimplify[trans,
  200.                      Dialogue -> sdialogue,
  201.                      Variables -> tlist],
  202.               trans ];
  203.         If [ InformUserQ[dialogue], Scan[explain, lt, Infinity] ];
  204.         lt ]
  205.  
  206. explain[ invlaplace[ f_, s_, rest__ ] ]:=
  207.     Message[ Transform::incomplete, "inverse Laplace transform", f, s ]
  208.  
  209. validvarQ[ x_Symbol ] := True
  210. validvarQ[ x_List ] := Apply[And, Map[VariableQ, x]]
  211. validvarQ[ x_ ] := False
  212.  
  213.  
  214. (*  M E S S A G E S  *)
  215.  
  216. InvLaPlace::badROC =
  217.     "Improper region of convergence in ``."
  218. InvLaPlace::badterms =
  219.     "Option Terms is not an integer (Terms -> ``)."
  220. InvLaPlace::novariables =
  221.     "No Laplace variables specified in inverse Laplace transform."
  222.  
  223.  
  224. (*  S T A T E     D E F I N I T I O N  *)
  225.  
  226. Lfactorfield = 1        (* factor denominator              *)
  227. Lexpandfield = 2        (* apply Expand[] to expression          *)
  228. Lexpandallfield = 3        (* apply ExpandAll[] to expression      *)
  229. Lthefunfield = 4        (* apply TheFunction to expression      *)
  230. Lapartfield = 5            (* apply Apart[] to expression          *)
  231. Lmyapartfield = 6        (* apply MyApart[] to rational poly.      *)
  232. Lnormalizefield = 7        (* normalize denominator          *)
  233. Lseriesfield = 8        (* approximate function with power series *)
  234. Lsplitfield = 9            (* split rat. fun. into rat. poly & other *)
  235. Ldefinitionfield = 10        (* apply defintion of inv. Laplace trans. *)
  236.  
  237. Lstatevariables = 10
  238.  
  239. initInvLstate[] := Table[True, {Lstatevariables}]
  240. nullInvLstate[] := Table[False, {Lstatevariables}]
  241.  
  242.  
  243. (*  S U P P O R T I N G     R O U T I N E S  *)
  244.  
  245. (*  Causal or anti-causal  *)
  246. leftsided[ p_, t_ ] :=
  247.     ( p /. { CStep[ b_. t + c_. ] -> CStep[ -b t + c ], t -> - t } )
  248.  
  249. rightsided[ p_, t_ ] := p
  250.  
  251. leftOrRightSided[ a_, p_, t_, rp_ ] := rightsided[ p, t ] /; InfinityQ[rp]
  252. leftOrRightSided[ a_, p_, t_, rp_ ] := leftsided[ p, t ]  /; SameQ[Re[a], rp]
  253. leftOrRightSided[ a_, p_, t_, rp_ ] := leftsided[ p, t ]  /; SameQ[-Re[a], rp]
  254. leftOrRightSided[ a_, p_, t_, rp_ ] := leftsided[ p, t ]  /; N[Re[a] > rp]
  255. leftOrRightSided[ a_, p_, t_, rp_ ] := rightsided[ p, t ]
  256.  
  257. (*  MyInvLaPlace  *)
  258. MyInvLaPlace[ e_, s_, t_, op_ ] :=
  259.     Block [ {rminus, rplus, valid},
  260.         rminus = SPSimplify[ GetRMinus[e] ];
  261.         rplus = SPSimplify[ GetRPlus[e] ];
  262.         valid = Apply[And, Thread[ToList[rminus] < ToList[rplus]]];
  263.         If [ TrueQ[! valid],
  264.              Message[ InvLaPlace::badROC, e ],
  265.              MyInvLaPlace[ TheFunction[e], s, t, rminus,
  266.                         rplus, initInvLstate[], op ] ] ] /;
  267.     LTransformQ[e]
  268.  
  269. MyInvLaPlace[ e_, s_, t_, op_ ] :=
  270.     MyInvLaPlace[ e, s, t, -Infinity, Infinity, initInvLstate[], op ]
  271.  
  272. MyInvLaPlace[ e_, s_, t_, rm_, rp_, st_, op_ ] :=
  273.     Block [    {laste, newe, newinvrules, norme, trace},
  274.  
  275.         (* add the time variable to the global list of time vars *)
  276.         AppendTo[tlist, t];
  277.  
  278.         (* generate the rules *)
  279.         newinvrules = TransformFixUp[ s, t, op, myinvlaplace, False,
  280.                           InvLaPlace, Null, Null ];
  281.         InvLaPlaceRules = Join[ newinvrules, OriginalInvLaPlaceRules ];
  282.  
  283.         (* determine the level of dialogue and pre-process f(t) *)
  284.         trace = SameQ[ Replace[Dialogue, op], All ];
  285.         norme = normalizeExpression[e, s];
  286.         If [ trace && ! SameQ[norme, e],
  287.              Print[e]; Print["is first rewritten as"];
  288.              Print[norme]; Print["Then, its inverse transform is"] ];
  289.  
  290.         (* take the 1-D inverse Laplace transform *)
  291.         newe = myinvlaplace[ norme, s, t, rm, rp, st, fixUp[op] ];
  292.         While [ ! SameQ[laste, newe],
  293.             If [ trace, Print[ newe ]; Print[ "which becomes" ] ];
  294.             laste = newe;
  295.             newe = MapAll[transform, laste] ];
  296.  
  297.         (* fix up the result *)
  298.         newe = postinvtransform[laste];
  299.         While [ ! SameQ[laste, newe],
  300.             If [ trace, Print[ newe]; Print[ "which becomes" ] ];
  301.             laste = newe;
  302.             newe = postinvtransform[laste] ];
  303.         If [ trace, Print[ newe] ];
  304.  
  305.         (* Simplification rule based on t being real variable *)
  306.         (* Mathematica 2.0 always assumes var to be complex   *)
  307.         If [ TrueQ[ $VersionNumber >= 2.0 ],
  308.              newe = newe /.  ( ((t^l_.) x_)^k_ :> t^(l k) x^k /;
  309.                     ( Positive[l] || Negative[l] ) &&
  310.                     ( Positive[k] || Negative[k] ) ) ];
  311.  
  312.         newe /. ( f_ Delta[t + a_.] :> (f /. t -> -a) Delta[t + a] ) ]
  313.  
  314.  
  315. (*  Kludge around the way Mathematica does Apart            *)
  316. (*  If the user has selected Apart -> All, then use approximate roots    *)
  317. (*  Otherwise, inverse transform it as a filter structure        *)
  318. (*  Enable denominator normalization just in case Apart doesn't do it    *)
  319.  
  320. breakUpRealValuedPoly[x_, s_, t_, rm_, rp_, st_, op_] :=
  321.     Block [    {dialogue, input, numer, state,},
  322.         state = SetStateField[st, Lmyapartfield, False];
  323.         state = SetStateField[state, Lnormalizefield, True];
  324.         dialogue = InformUserQ[op];
  325.         If [ dialogue,
  326.              Print[ "( since the Apart option is not All," ];
  327.              Print[ "  the inverse Laplace transform of" ];
  328.              Print[ "  ", x ];
  329.              Print[ "  will be an IIR filter whose input is" ];
  330.              Print[ "  the inverse transform of the numerator. )" ] ];
  331.  
  332.         numer = Numerator[x];
  333.         input = If [ realValuedPolynomialQ[numer, s] &&
  334.                      realValuedCoefficientsQ[numer, s],
  335.                  CFIR[t, CoefficientList[numer, s]],
  336.                  myinvlaplace[numer, s, t, rm, rp, state, op] ];
  337.         CIIR[t, CoefficientList[Denominator[x], s]][input] ] /;
  338.     ! SameQ[ Replace[Apart, op], All ]
  339.  
  340. breakUpRealValuedPoly[x_, s_, t_, rm_, rp_, st_, op_] :=
  341.     partialFractionsKludge[x, s, t, rm, rp, st, op, N, "approximate"]
  342.  
  343. (*  definitionDialogue  *)
  344. definitionDialogue[ oldexpr_, trans_, operator_, options_ ] :=
  345.     Block [    {},
  346.         If [ InformUserQ[options],
  347.              Print[ "( after using ", operator,
  348.                 "  to apply the definition to" ];
  349.              Print[ "  ", oldexpr ];
  350.              Print[ "  to find its transform"];
  351.              Print[ "  ", trans, " )" ] ];
  352.         trans ]
  353.  
  354. (*  fixUp  *)
  355. fixUp[ op_ ] := { Apart -> Replace[Apart, op],
  356.           Definition -> Replace[Definition, op],
  357.           Dialogue -> Replace[Dialogue, op],
  358.           Simplify -> Replace[Simplify, op],
  359.           Terms -> Replace[Terms, op] }
  360.  
  361.  
  362. (*  multiDInvLaPlace  *)
  363. multiDInvLaPlace[ e_, s_, t_, op_ ] :=
  364.     MultiDInvTransform [ e, s, t, op, LTransformQ,
  365.                  MyInvLaPlace, MakeLObject, Global`t ]
  366.  
  367. (*  normalizePoly and normalizeExpression  *)
  368. truepoly[x_, s_] := (! FreeQ[x, s]) && PolynomialQ[x, s]
  369.  
  370. unnormalizedq[x_, s_] := ! SameQ[1, Coefficient[x, s, Exponent[x, s]]]
  371.  
  372. normalizePoly[poly_, s_, k_] :=
  373.     Block [    {coeff, maxpower},
  374.         maxpower = Exponent[poly, s];
  375.         coeff = Coefficient[poly, s, maxpower];
  376.         coeff^k ( Distribute[poly / coeff]^k ) ]
  377.  
  378. normalizeExpression[x_, s_] :=
  379.     Block [ {},
  380.         x /. { (poly_Plus)^k_. :> normalizePoly[Expand[poly], s, k] /;
  381.              truepoly[poly, s],
  382.                poly_ :> Expand[poly] /; truepoly[poly, s] } ]
  383.  
  384. (*  normalizeRatPoly  *)
  385. normalizeRatPoly[x_, s_] := normalizeRatPoly[Numerator[x], Denominator[x], s]
  386.  
  387. normalizeRatPoly[numer_, denom_^k_., s_] :=
  388.     Block [ {coeff, maxpower},
  389.         maxpower = Exponent[denom, s];
  390.         coeff = Coefficient[denom, s, maxpower];
  391.         coeff^k numer / ( Distribute[denom / coeff]^k ) ]
  392.  
  393. (*  partialFractionsDialogue  *)
  394. partialFractionsDialogue[p_, s_, options_] :=
  395.     partialFractionsDialogue[ p, s, options,
  396.                   normalizeExpression[Apart[p, s], s] ]
  397.  
  398. partialFractionsDialogue[p_, s_, options_, partfrac_, rootinfo_:"exact"] :=
  399.     Block [ {dialogue},
  400.         dialogue = InformUserQ[options];
  401.  
  402.         (*  Tell the user that we're applying partial fractions  *)
  403.         If [ dialogue && ! SameQ[p, partfrac],
  404.              Print["( after performing partial fraction expansion on"];
  405.              Print["  ", p];
  406.              Print["  using its ", rootinfo, " roots:" ];
  407.              Print["  ", partfrac, " . )" ] ];
  408.  
  409.         partfrac ]
  410.  
  411. (*  partialFractionsKludge -- use MyApart because Apart has failed us  *)
  412. partialFractionsKludge[ x_, s_, t_, rm_, rp_, st_, op_,
  413.             filter_:Identity, roots_:"exact" ] :=
  414.     Block [    {partfrac, result, state},
  415.         state = SetStateField[st, Lapartfield, False];
  416.         state = SetStateField[state, Lmyapartfield, False];
  417.         state = SetStateField[state, Lnormalizefield, False];
  418.         partfrac = normalizeExpression[MyApart[x, s, filter], s];
  419.         result = partialFractionsDialogue[x, s, op, partfrac, roots];
  420.         myinvlaplace[ result, s, t, rm, rp, state, op ] ]
  421.  
  422. (*  postinvtransform  *)
  423. postinvtransform[f_] := ( f /. postinvrules )
  424.  
  425. postinvrules = {
  426.     derivative[f_, t_Symbol, k_, False] :>
  427.         D[ postinvtransform[f], {t, k} ],
  428.     derivative[f_, t_Symbol, k_, v_] :>
  429.         Block [    {curderivative, ffun, i, limit, result},
  430.             ffun = postinvtransform[f];
  431.             result = Unit[k-1][t] Limit[ffun, t -> v];
  432.             curderivative = ffun;
  433.             For [ i = 2, i <= k, i++,
  434.                   curderivative = D[curderivative, t];
  435.                   limit = Limit[curderivative, t -> v];
  436.                   result += Unit[k - i][t] limit ];
  437.             result + D[curderivative, t] ],
  438.     integrate[a_ + b_, t_, dummy_Symbol] :>
  439.         integrate[postinvtransform[a], t, dummy] +
  440.         integrate[postinvtransform[b], t, dummy],
  441.     integrate[f_. CStep[dummy_ + b_.], t_, dummy_Symbol ] :>
  442.         Integrate[f, {dummy, -b, t}] CStep[t + b],
  443.     integrate[f_. CStep[a_. dummy_ + b_.], t_, dummy_Symbol ] :>
  444.         ( CStep[a t + b] / a ) *
  445.         Integrate[f /. t -> t / a,
  446.               {dummy, Max[-b, -Infinity a], a t}],
  447.     substitutefort[f_, t_, newt_] :>
  448.         ( postinvtransform[f] /. t -> newt ),
  449.  
  450.     (*  These two rules are needed to support a series     *)
  451.     (*  expansion which returns a constant term plus    *)
  452.     (*  terms of the form  approx * s^r   .            *)
  453.     myinvlaplace[ a_, s_, t_, rm_, rp_, st_, op_] :> a Delta[t] /;
  454.         FreeQ[a, s],           (* Converges all ROC *)
  455.  
  456.     myinvlaplace[ a_. s_^r_., s_, t_, rm_, rp_, st_, op_] :>
  457.         a Unit[r][t] /;
  458.         FreeQ[{a,r}, s],       (* Converges all ROC *)
  459.  
  460.     (*  Attempt a series expansion about s = 0        *)
  461.     (*  This is the strategy when all else has failed.    *)
  462.     (*  Please keep this as the last rule.            *)
  463.     myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
  464.         SeriesStrategy[ x, s, t, rm, rp, st, op ] /;
  465.         GetStateField[st, Lseriesfield],
  466.  
  467.     x_ :> x
  468. }
  469.  
  470. (*  realAndPosPart  *)
  471. realposcheckq[x_] := Implies[ NumberQ[x], RealQ[x] && x > 0 ]
  472. realnegcheckq[x_] := Implies[ NumberQ[x], RealQ[x] && x < 0 ]
  473.  
  474. realAndPosPart[0] := 0
  475. realAndPosPart[x_Times] := Map[realAndPosPart, x]
  476. realAndPosPart[Complex[a_Integer, b_Integer]] := realAndPosPart[GCD[a,b]]
  477. realAndPosPart[Complex[b_, b_]] := realAndPosPart[b]
  478. realAndPosPart[Complex[0, b_]] := realAndPosPart[b]
  479. realAndPosPart[x_Symbol] := x
  480. realAndPosPart[x_] := x /; realposcheckq[ N[x] ]
  481. realAndPosPart[x_] := -x /; realnegcheckq[ N[x] ]
  482. realAndPosPart[x_] := 1
  483.  
  484. (*  realValuedCoefficientsQ  *)
  485. realValuedCoefficientsQ[z_, t_] :=
  486.     Apply[Or, Map[RealQ, CoefficientList[z, t]]]
  487.  
  488. (*  realValuedPolynomialQ  *)
  489. realValuedPolynomialQ[z_] := realValuedPolynomialQ[z, Global`x]
  490.  
  491. realValuedPolynomialQ[c_. t_^r_., t_] := IntegerQ[r] && r > 0 && RealValuedQ[c]
  492. realValuedPolynomialQ[c_, t_] := RealValuedQ[c]
  493. realValuedPolynomialQ[x_ + y_, t_] :=
  494.     realValuedPolynomialQ[x, t] && realValuedPolynomialQ[y, t]
  495.  
  496. (*  seriesDialogue  *)
  497. seriesDialogue[x_, expansion_, terms_, options_] :=
  498.     Block [    {dialogue},
  499.         dialogue = InformUserQ[options];
  500.  
  501.         If [ dialogue,
  502.              Print["( after breaking up the expression"];
  503.              Print["  ", x];
  504.              Print["  into its series representation"];
  505.              Print["  ", expansion];
  506.              Print["  the inverse Laplace transform will now be",
  507.                "  applied to the first ", terms, " terms. )"] ] ]
  508.  
  509. (*  SeriesExpansion  *)
  510. SeriesExpansion[ True, x_, s_, power_ ] := Series[x, {s, 0, power}]
  511.  
  512. SeriesExpansion[ False, x_, s_, power_ ] :=
  513.     Block [    {negx, expansion},
  514.         negx = x /. s -> s^-1;
  515.         expansion = Series[negx, {s, 0, power}];
  516.         expansion /. s -> s^-1 ]
  517.  
  518. (*  SeriesStrategy  *)
  519. SeriesStrategy[ x_, s_, t_, rm_, rp_, st_, op_ ] :=
  520.     Block [ {expansion, exponents, posexpandflag, state, terms},
  521.         terms = Replace[Terms, op];
  522.         state = SetStateField[st, Lseriesfield, False];
  523.  
  524.         (*  If the Terms option is disabled, then take no action.  *)
  525.         If [ ! terms,
  526.              Return[ myinvlaplace[x, s, t, rm, rp, state, op] ] ];
  527.  
  528.         If [ ! IntegerQ[terms],
  529.              Message[ InvLaPlace::badterms, terms ];
  530.              Return[ myinvlaplace[x, s, t, rm, rp, state, op] ] ];
  531.  
  532.         (*  See if we should expand in pos. or neg. powers    *)
  533.         exponents = GetAllExponents[x, s];
  534.         posexpandflag = If [ Apply[And, Thread[exponents >= 0]],
  535.                      True, False, False ];
  536.  
  537.         (*  First expansion  *)
  538.         expansion = SeriesExpansion[posexpandflag, x, s, terms - 1];
  539.  
  540.         (*  Second expansion if first failed  *)
  541.         If [ ! SameQ[Head[expansion], SeriesData],
  542.              expansion = SeriesExpansion[Not[posexpandflag], x, s, terms - 1] ];
  543.  
  544.         (*  If this fails, call MyInvLaPlace again        *)
  545.         (*  but with series expansion disabled            *)
  546.         If [ ! SameQ[Head[expansion], SeriesData],
  547.              myinvlaplace[ x, s, t, rm, rp, state, op ],
  548.              seriesDialogue[x, expansion, terms, op];
  549.              state = nullInvLstate[];
  550.              Map [ myinvlaplace[#1, s, t, rm, rp, state, op]&, 
  551.                Normal[expansion] ] ] ]
  552.  
  553. (*  similarityQ  *)
  554. similarityScale[f_, s_] := realAndPosPart[ ScalingFactor[f, s] ]
  555. similarityQ[f_, s_] :=
  556.     Block [    {scale},
  557.         scale = similarityScale[f, s];
  558.         (! SameQ[scale, 1]) && (! SameQ[scale, 0]) ]
  559.  
  560. (*  splitratfun *)
  561. keepterm[ x_ + y_, s_ ] := keepterm[x, s] + keepterm[y, s]
  562. keepterm[ a_. s_^n_., s_ ] := a s^n /; FreeQ[{a, n}, s]
  563. keepterm[ a_, s_ ] := a /; FreeQ[a, s]
  564. keepterm[ a_, s_ ] := 0
  565.  
  566. splitratfun[x_, s_] := splitratfun[ Numerator[x], Denominator[x], s ]
  567.  
  568. splitratfun[ numer_, denom_, s_ ] :=
  569.     Block [ {expnumer, polyexpr, restexpr},
  570.         expnumer = Expand[numer];
  571.         polyexpr = keepterm[expnumer, s];
  572.         restexpr = expnumer - polyexpr;
  573.         { polyexpr / denom, restexpr / denom } ]
  574.  
  575. (*  transform  *)
  576. transform[ expr_ ] :=
  577.     If [ SameQ[Head[expr], myinvlaplace],
  578.          Replace[expr, InvLaPlaceRules],
  579.          expr ]
  580.  
  581.  
  582. Format[ myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] ] :=
  583.     SequenceForm[ ColumnForm[{"L" Superscript[-1],
  584.                   "  " ~StringJoin~ ToString[s]}],
  585.               { x } ]
  586.  
  587. Format[ derivative[x_, t_Symbol, k_, flag_] ] := 
  588.     SequenceForm[ ColumnForm[{"D" Superscript[k],
  589.                   "  " ~StringJoin~ ToString[t]}],
  590.               { x } ]
  591.  
  592. Format[    substitutefort[f_, t_, newt_] ] :=
  593.     SequenceForm[ {f}, Subscript[t -> newt] ]
  594.  
  595. Format[ integrate[f_. CStep[a_. dummy_ + b_.], t_, dummy_Symbol ] ] :=
  596.     SequenceForm[ "Integrate[", f, ", ", {dummy, 0, t}, "] ",
  597.               CStep[a t + b] ]
  598.  
  599. Format[ integrate[f_, t_, dummy_Symbol ] ] :=
  600.     SequenceForm[ "Integrate[", f, ", ", {dummy, 0, t}, "]" ]
  601.  
  602.  
  603. (*  B E G I N     R U L E S  *)
  604.  
  605.  
  606. InvLaPlaceRules = {}
  607.  
  608.  
  609. OriginalInvLaPlaceRules = { 
  610.  
  611.  
  612.   (* I.   Rational Transform Pairs                *)
  613.  
  614.   (*      A.  constant functions                *)
  615.   myinvlaplace[ a_, s_, t_, rm_, rp_, st_, op_] :> a Delta[t] /;
  616.     FreeQ[a, s],                    (* Converges all ROC *)
  617.  
  618.  
  619.   (*      B.  constant times exponential            *)
  620.   myinvlaplace[ a_. Exp[b_. s_], s_, t_, rm_, rp_, st_, op_ ] :>
  621.     a Delta[t + b] /;
  622.     FreeQ[{a, b}, s],
  623.  
  624.  
  625.   (*      C.  inverse of sinusoidal forms            *)
  626.   myinvlaplace[ c_. ( w_. Cos[w_] + s_ Sin[w_] ) / (s_^2 + w_^2), s_, t_,
  627.         rm_, rp_, st_, op_ ] :>
  628.     leftOrRightSided[ 0, c Sin[w + w t] CStep[t], t, rp ] /;
  629.     FreeQ[{c, w}, s],
  630.  
  631.   myinvlaplace[ c_. ( s_ Cos[b_] - w_. Sin[b_] ) / ( s_^2 + w_^2 ), s_, t_,
  632.         rm_, rp_, st_, op_ ] :>
  633.     leftOrRightSided[ 0, c Cos[b + w t] CStep[t], t, rp ] /;
  634.     FreeQ[{w,b,c}, s],
  635.  
  636.  
  637.   (*      D.  sections that are functions of single power of s    *)
  638.   (*          1.  single pole sections                *)
  639.   myinvlaplace[ b_. (s_ + a_.)^n_., s_, t_, rm_, rp_, st_, op_ ] :>
  640.     leftOrRightSided[ a,
  641.               b t^-(n+1) Exp[- a t] CStep[t] / Gamma[-n],
  642.               t, rp ] /;
  643.     FreeQ[{a,b,n}, s] && Negative[n],
  644.  
  645.   myinvlaplace[ b_. / (s_ + a_.)^n_., s_, t_, rm_, rp_, st_, op_ ] :>
  646.     leftOrRightSided[ a,
  647.               b t^(n-1) Exp[- a t] CStep[t] / Gamma[n],
  648.               t, rp ] /;
  649.     FreeQ[{a,b,n}, s],
  650.  
  651.   (*          2.  single zero sections                *)
  652.   myinvlaplace[ b_. (s_ + a_.)^n_., s_, t_, rm_, rp_, st_, op_ ] :>
  653.     b Exp[- a t] Unit[n][t] /;
  654.     FreeQ[{a,b,n}, s],
  655.  
  656.  
  657.  
  658.   (*      D.  second-order sections                *)
  659.   (*          only encode the cosine/sine pairs because Mma    *)
  660.   (*          will automatically convert Sin[I x] to I Sinh[x]    *)
  661.   (*          and Cos[I x] to Cosh[x].                *)
  662.   myinvlaplace[ b_. / ( s_^2 + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
  663.     leftOrRightSided[Sqrt[a], b Sin[Sqrt[a] t] CStep[t] / Sqrt[a], t, rp] /;
  664.     FreeQ[{a,b}, s],
  665.  
  666.   myinvlaplace[ b_. s_ / ( s_^2 + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
  667.     leftOrRightSided[ Sqrt[a], b Cos[Sqrt[a] t] CStep[t], t, rp ] /; 
  668.     FreeQ[{a,b}, s],
  669.  
  670.   myinvlaplace[ (s_ + a_)/(s_^2 + c_. s_ + d_), s_, t_, rm_, rp_, st_, op_ ] :>
  671.     Block [    {b, invtrans},
  672.         b = Sqrt[d - a^2];            (* d = a^2 + b^2 *)
  673.         invtrans = Exp[-a t] Cos[b t] CStep[t];
  674.         leftOrRightSided[-a, invtrans, t, rp] ] /;
  675.     FreeQ[{a,c,d}, s] && SameQ[c, 2 a],
  676.  
  677.   myinvlaplace[ k_. / ( s_^2 + c_. s_ + d_ ), s_, t_, rm_, rp_, st_, op_ ] :>
  678.     Block [    {a},
  679.         a = c/2;        (* complete the square *)
  680.         Exp[-a t] myinvlaplace[ k / ( s^2 + d - a^2 ), s, t,
  681.                     rm, rp, st, op ] ] /;
  682.     FreeQ[{c,d,k}, s],
  683.  
  684.  
  685.   (*      E.  third order sections                *)
  686.   myinvlaplace[ b_. /( s_^3 + a_), s_, t_, rm_, rp_, st_, op_ ] :>
  687.     leftOrRightSided[ Max[-Re[(-a)^(1/3)/2], -Re[a^(1/3)]],
  688.               b / (3 a^(2/3)) (Exp[- a^(1/3) t] - Exp[a^(1/3) t/2]
  689.                 ( Cos[ a^(1/3) t Sqrt[3] / 2 ] - 
  690.                   Sqrt[3] Sin[ a^(1/3) t Sqrt[3] / 2 ] ) )
  691.               CStep[t],
  692.               t,
  693.               rp ] /;
  694.     FreeQ[{a, b}, s],                (* [Churchill] *)
  695.  
  696.  
  697.   (*      F.  fourth order sections                *)
  698.   myinvlaplace[ b_. / ( s_^2 + a_ )^2, s_, t_, rm_, rp_, st_, op_ ] :>
  699.     Block [    {k},
  700.         k = Sqrt[a] t;
  701.         leftOrRightSided[ Sqrt[-a],
  702.                      b (Sin[k] - k Cos[k]) CStep[t] / (2 a^(3/2)),
  703.                   t, rp ] ] /;
  704.     FreeQ[{a, b}, s],
  705.  
  706.   myinvlaplace[ b_. s_ / ( s_^2 + a_ )^2, s_, t_, rm_, rp_, st_, op_ ] :>
  707.     leftOrRightSided[ Sqrt[-a],
  708.               b t Sin[Sqrt[a] t] CStep[t] / (2 Sqrt[a]),
  709.               t, rp ] /;
  710.     FreeQ[{a, b}, s],
  711.  
  712.   myinvlaplace[ b_. / ( s_^4 + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
  713.     Block [    {k},
  714.         k = (a/4)^(1/4);
  715.         leftOrRightSided[ Re[ a^(1/4) / Sqrt[2] ],
  716.                   b (Sin[k t] Cosh[k t] - Cos[k t] Sinh[k t])
  717.                     CStep[t] / ( 4 k^3 ),
  718.                   t,
  719.                   rp ] ] /;
  720.     FreeQ[{a, b}, s],
  721.  
  722.   myinvlaplace[ b_. s_ / ( s_^4 + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
  723.     Block [    {k},
  724.         k = (a/4)^(1/4);
  725.         leftOrRightSided[ Re[ a^(1/4) / Sqrt[2] ], 
  726.                   b Sin[k t] Sinh[k t] CStep[t] / ( 2 k^2 ),
  727.                   t,
  728.                   rp ] ] /;
  729.      FreeQ[{a, b}, s],
  730.  
  731.  
  732.   (*      G.  higher order sections                *)
  733.   myinvlaplace[ b_. s_^n_., s_, t_, rm_, rp_, st_, op_ ] :>
  734.       leftOrRightSided[ 0, b t^-(n + 1) / Gamma[-n] CStep[t], t, rp ] /; 
  735.     Negative[n] && FreeQ[b, s],
  736.  
  737.   myinvlaplace[ b_. s_^n_, s_, t_, rm_, rp_, st_, op_ ] :>
  738.     Block [    {index = Unique["index"]}, 
  739.         leftOrRightSided[ 0,
  740.                b CStep[t] 2^(-n-1/2) t^(-n-1) /
  741.                ( Product[index, {index, 1, -2 n - 1}] Sqrt[Pi] ),
  742.                t, rp ] ] /; 
  743.     IntegerQ[ - n - 1/2 ] && FreeQ[b, s],
  744.  
  745.   myinvlaplace[ b_. ( s_^2 + a_ )^k_., s_, t_, rm_, rp_, st_, op_ ] :>
  746.     leftOrRightSided[ a,
  747.               b Sqrt[Pi] t^(-k - 1/2) / ((2 Sqrt[a])^(-k - 1/2))
  748.                 BesselJ[-k - 1/2, Sqrt[a] t ] CStep[t] / Gamma[-k],
  749.               t, rp ] /;
  750.     FreeQ[{a, b}, s] && Negative[k],
  751.  
  752.   myinvlaplace[ b_. / ( s_^2 + a_ )^k_., s_, t_, rm_, rp_, st_, op_ ] :>
  753.     leftOrRightSided[ a,
  754.               b Sqrt[Pi] t^(k -1/2) / ((2 Sqrt[a])^(k - 1/2))
  755.                 BesselJ[k - 1/2, Sqrt[a] t ] CStep[t] / Gamma[k],
  756.               t, rp ] /;
  757.     FreeQ[{a,b,k}, s],
  758.  
  759.  
  760.   (*      H.  partial fractions of rational polynomials        *)
  761.  
  762.   (*          1. split expression into rational polys and rest  *)
  763.   (*             this must occur before partial fractions    *)
  764.   (*         because Apart ONLY handles ratio of polys    *)
  765.   myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
  766.     Block [    {splitlist, state},
  767.         splitlist = splitratfun[x, s];
  768.         state = SetStateField[st, Lsplitfield, False];
  769.         myinvlaplace[ First[splitlist], s, t, rm, rp, state, op] +
  770.           myinvlaplace[ Second[splitlist], s, t, rm, rp, state, op] ] /;
  771.     GetStateField[st, Lsplitfield] && GetStateField[st, Lapartfield] &&
  772.       truepoly[Denominator[x], s] && (! PolynomialQ[Numerator[x], s]),
  773.  
  774.   (*          2. Reduce all rational polynomials into        *)
  775.   (*             a sum of single pole, second order sections,    *)
  776.   (*             and sections like 1 / ( s^4 + a^4 )        *)
  777.   (*             Enable denominator normalization just in case    *)
  778.   (*             partial fractions does not do it automatically    *)
  779.   myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
  780.     Block [    {state},
  781.         state = SetStateField[st, Lapartfield, False];
  782.         state = SetStateField[state, Lnormalizefield, True];
  783.         myinvlaplace[ partialFractionsDialogue[x, s, op], s, t,
  784.                   rm, rp, state, op ] ] /;
  785.     GetStateField[st, Lapartfield] && RationalPolynomialQ[x, s],
  786.  
  787.  
  788.   (*      I.  partial fractions that Apart could not resolve    *)
  789.   (*          perform MyApart to redo partial fractions for    *)
  790.   (*          denominators with real-valued coefficients    *)
  791.   (*          that are not handled by Apart.         *)
  792.   myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
  793.     breakUpRealValuedPoly[x, s, t, rm, rp, st, op] /;
  794.     GetStateField[st, Lmyapartfield] &&
  795.       realValuedPolynomialQ[Denominator[x], s] &&
  796.       realValuedCoefficientsQ[Denominator[x], s],
  797.  
  798.  
  799.  
  800.  
  801.   (* II.  Non-rational Transform Pairs                *)
  802.  
  803.   (*      A.  square root forms                    *)
  804.   myinvlaplace[ 1 / Sqrt[s_ + b_.], s_, t_, rm_, rp_, st_, op_ ] :>
  805.     leftOrRightSided[ 0, Exp[- b t] CStep[t] / (Sqrt[t] Sqrt[Pi]),
  806.               t, Infinity ] /;
  807.     FreeQ[b, s],
  808.  
  809.   myinvlaplace[ s_ / ( s_ + a_ )^(3/2), s_, t_, rm_, rp_, st_, op_ ] :>
  810.     leftOrRightSided[ -Re[a], Exp[- a t] (1 - 2 a t) CStep[t] / Sqrt[Pi],
  811.               t, rp ] /;
  812.     FreeQ[a, s],
  813.  
  814.   myinvlaplace[ 1 / ((s_ + a_.)^(3/2) Sqrt[s_ + b_.]), s_, t_,
  815.         rm_, rp_, st_, op_ ] :>
  816.     leftOrRightSided[ a,
  817.         t Exp[- (a + b) t / 2] ( BesselI[0, (a - b) t / 2] -
  818.                      BesselI[1, (a - b) t / 2] ) CStep[t],
  819.         t, rp ] /;
  820.     FreeQ[{a,b}, s] && ! SameQ[a, b],
  821.  
  822.   myinvlaplace[ 1 / ( Sqrt[s_] + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
  823.     (1 / (Sqrt[Pi] Sqrt[t]) - a Exp[a^2 t](1 - Erf[a Sqrt[t]])) CStep[t] /; 
  824.     FreeQ[a, s],
  825.  
  826.   (*          general form of #38 and #39 [Churchill, 326],    *)
  827.   (*          Sqrt[s] / ( s + a ), augmented by s -> s + b    *)
  828.   myinvlaplace[ Sqrt[s_ + b_.] / ( s_ + c_ ), s_, t_, rm_, rp_, st_, op_ ] :>
  829.       leftOrRightSided[c,
  830.              ( 1 / (Sqrt[Pi] Sqrt[t]) - Sqrt[b - c] Exp[(b - c) t]
  831.                Erf[ Sqrt[(b - c) t] ] ) CStep[t],
  832.              t, rp ] /;
  833.     FreeQ[{b,c}, s],
  834.  
  835.   (*          general form of #40 and #41 [Churchill, 326],    *)
  836.   (*          1 / ( Sqrt[s] (s + a) ), augmented by s -> s + b    *)
  837.   myinvlaplace[ 1 / (Sqrt[s_ + b_.] (s_ + c_)), s_, t_, rm_, rp_, st_, op_ ] :>
  838.     leftOrRightSided[c,
  839.              Exp[- c t] Erf[ Sqrt[(b - c) t] ] CStep[t] /
  840.                Sqrt[b - c],
  841.              t, rp ] /;
  842.     FreeQ[{b,c}, s],
  843.  
  844.   (*          forms that inverse transform to Bessel functions    *)
  845.   (*          1/Sqrt[s^2 + a] covered as a higher order section    *)
  846.   myinvlaplace[ ( Sqrt[s_^2 + a_] - s_ )^k_, s_, t_, rm_, rp_, st_, op_ ] :>
  847.     leftOrRightSided[ a, k Sqrt[a]^k BesselJ[k, Sqrt[a] t] CStep[t] / t,
  848.               t, rp ] /;
  849.     FreeQ[{a,k}, s] && Positive[k],
  850.  
  851.   myinvlaplace[ ( Sqrt[s_^2 + a_] - s_ )^v_. / Sqrt[s_^2 + a_], s_, t_,
  852.         rm_, rp_, st_, op_ ] :>
  853.     Block [    {timefun},
  854.         Assuming[ Positive[v + 1], op ];
  855.         timefun = Sqrt[a]^v BesselJ[v, Sqrt[a] t] CStep[t];
  856.         leftOrRightSided[a, timefun, t, rp] ] /;
  857.     FreeQ[{a,v}, s] && Implies[NumberQ[N[v]], N[v > -1]],
  858.  
  859.   myinvlaplace[ 1 / ( Sqrt[s_ + a_.] Sqrt[s_ + b_.] ), s_, t_,
  860.         rm_, rp_, st_, op_ ] :>
  861.     Exp[- (a + b) t / 2] BesselI[0, (a - b) t / 2] CStep[t] /;
  862.     FreeQ[{a, b}, s],
  863.  
  864.   myinvlaplace[ ( s_ - Sqrt[s_^2 + a_] )^v_. / Sqrt[s_^2 + a_], s_, t_,
  865.         rm_, rp_, st_, op_ ] :>
  866.     Block [    {timefun},
  867.         Assuming[ Positive[v + 1], op];
  868.         timefun = Sqrt[-a]^v BesselI[v, Sqrt[-a] t] CStep[t];
  869.         leftOrRightSided[ a, timefun, t, rp ] ] /;
  870.     FreeQ[{a,v}, s] && Implies[NumberQ[N[v]], N[v > -1]],
  871.  
  872.  
  873.   (*      B.  hyperbolic forms                     *)
  874.   myinvlaplace[ Tanh[k_. s_] / s_, s_, t_, rm_, rp_, st_, op_ ] :>
  875.     Block [    {n, context}, 
  876.         context = $Context;
  877.         $Context = "Global`";
  878.         n = Unique["n"];            (* n is an integer *)
  879.         $Context = context;
  880.         (-1)^(n-1) ( CStep[t - 2 k (n - 1)] - CStep[t - 2 k n] ) ] /;
  881.     FreeQ[k, s],
  882.  
  883.   myinvlaplace[ Coth[b_. s_] / ( s_^2 + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
  884.     Abs[Sin[Sqrt[a] t] ] / Sqrt[a] CStep[t] /;
  885.     b == ( Pi / ( 2 Sqrt[a] ) ),
  886.  
  887.  
  888.   (*      C.  exponential forms                    *)
  889.   myinvlaplace[ Exp[a_. Sqrt[s_]] / Sqrt[s_], s_, t_,
  890.         rm_, rp_, st_, op_ ] :>
  891.     Exp[-a^2 / (4 t)] CStep[t] / (Sqrt[Pi] Sqrt[t]) /;
  892.     FreeQ[a, s],
  893.  
  894.   myinvlaplace[ Exp[k_. / s_] s_^u_, s_, t_, rm_, rp_, st_, op_ ] :>
  895.     Block [ {arg, defarg, exp, order},
  896.         defarg = 2 Sqrt[k t];
  897.         arg = 2 Sqrt[k] Sqrt[t];
  898.         order = -u - 1;
  899.         exp = order/2;
  900.         Which [ Negative[k],
  901.               t^exp BesselJ[order, arg] CStep[t] / k^exp,
  902.             Positive[k],
  903.               t^exp BesselI[order, arg] CStep[t] / k^exp,
  904.             True,
  905.               t^exp BesselI[order, defarg] CStep[t] / k^exp ] ] /;
  906.     Negative[u],
  907.  
  908.   myinvlaplace[ Exp[k_. / s_], s_, t_, rm_, rp_, st_, op_ ] :>
  909.     Block [    {arg, scale},
  910.         arg = 2 Sqrt[-k] Sqrt[t];
  911.         scale = Sqrt[-k] / ( 2 Sqrt[t] );
  912.         Delta[t] + scale ( BesselJ[-1, arg] - BesselJ[1, arg] ) ] /;
  913.     FreeQ[k, s],
  914.  
  915.   myinvlaplace[ Exp[ k_ Sqrt[s_] ], s_, t_, rm_, rp_, st_, op_ ] :>
  916.     - k Exp[- k^2 / (4 t)] CStep[t] / (2 Sqrt[Pi] Sqrt[t^3]) /;
  917.     Negative[k],
  918.  
  919.  
  920.   (*      D.  logarithmic forms                    *)
  921.   myinvlaplace[ Log[c_. s_ + b_.] (c_. s_ + b_.)^k_., s_, t_,
  922.         rm_, rp_, st_, op_ ] :>
  923.        Exp[- b t / Abs[c]] (t / Abs[c])^(-k - 1)
  924.       ( Gamma[-k] PolyGamma[-k] / Abs[Gamma[-k]]^2 -
  925.         Log[t / Abs[c]] / Gamma[-k] ) CStep[t] / Abs[c] /;
  926.     FreeQ[{b,c}, s] && Negative[k],
  927.  
  928.   myinvlaplace[ Log[c_. s_ + b_.] / (c_. s_ + b_.)^k_., s_, t_,
  929.         rm_, rp_, st_, op_ ] :>
  930.        Exp[- b t / Abs[c]] (t / Abs[c])^(k - 1)
  931.       ( Gamma[k] PolyGamma[k] / Abs[Gamma[k]]^2 -
  932.         Log[t / Abs[c]] / Gamma[k] ) CStep[t] / Abs[c] /;
  933.     FreeQ[{b,c,k}, s],
  934.  
  935.   myinvlaplace[ Log[b_. s_ + a_] / s_, s_, t_, rm_, rp_, st_, op_ ] :>
  936.     ( Log[a] - ExpIntegralEi[- a t / Abs[b]] ) CStep[t] /;
  937.     FreeQ[{a,b}, s],
  938.  
  939.   myinvlaplace[ Log[b_. s_ + a_.], s_, t_, rm_, rp_, st_, op_ ] :>
  940.       - Abs[b] Exp[- a t / Abs[b]] CStep[t] / t /;
  941.     FreeQ[{a,b}, s],
  942.  
  943.   myinvlaplace[ Log[a_. + b_. s_^2 ] / s_, s_, t_, rm_, rp_, st_, op_ ] :>
  944.     2 ( Log[Sqrt[a]] - CosIntegral[Sqrt[a] t / Sqrt[b]] ) CStep[t] /;
  945.      FreeQ[{a,b}, s],
  946.  
  947.   myinvlaplace[ Log[b_. s_] / (c_. s_^2 + 1), s_, t_, rm_, rp_, st_, op_ ] :>
  948.     Block [    {time = t / Abs[b]},
  949.         ( Cos[time] SinIntegral[time] -
  950.           Sin[time] CosIntegral[time] ) CStep[t] / Abs[b] ] /;
  951.     FreeQ[{b,c}, s] && c == b^2,
  952.  
  953.   myinvlaplace[ s_ Log[b_. s_] / (c_. s_^2 + 1), s_, t_, rm_, rp_, st_, op_ ] :>
  954.     Block [    {time = t / Abs[b]},
  955.         ( - Cos[time] CosIntegral[time] -
  956.             Sin[time] SinIntegral[time] ) CStep[t] / (b Abs[b]) ] /;
  957.     FreeQ[{b,c}, s] && c == b^2,
  958.  
  959.   myinvlaplace[ Log[ a_ / b_ ], s_, t_, rm_, rp_, st_, op_ ] :>
  960.       myinvlaplace[ Log[a] - Log[b], s, t, rm, rp, st, op ],
  961.  
  962.   myinvlaplace[ Log[ a_ b_ ], s_, t_, rm_, rp_, st_, op_ ] :>
  963.       myinvlaplace[ Log[a] + Log[b], s, t, rm, rp, st, op ],
  964.  
  965.  
  966.   (*      E.  arc tangent                    *)
  967.   myinvlaplace[ ArcTan[ k_ / s_ ], s_, t_, rm_, rp_, st_, op_ ] :>
  968.       Sin[ k t ] CStep[t] / t /; FreeQ[k,s],
  969.  
  970.  
  971.   (*      F.  BesselK                        *)
  972.   myinvlaplace[ BesselK[ 0, k_. s_ ], s_, t_, rm_, rp_, st_, op_ ] :>
  973.        1 / Sqrt[t^2 - k^2] CStep[t - k] /;
  974.     FreeQ[k, s],
  975.  
  976.   myinvlaplace[ BesselK[-1, v_. Sqrt[s_]] / Sqrt[s_], s_, t_,
  977.         rm_, rp_, st_, op_ ] :>
  978.     Exp[-v^2/(4 t)] CStep[t] / v /;
  979.     FreeQ[v, s],
  980.  
  981.  
  982.   (* III. Transform Properties                    *)
  983.  
  984.  
  985.   (*      A.  homogeneity -- pick off constants            *)
  986.   myinvlaplace[ a_ f_, s_, t_, rm_, rp_, st_, op_ ] :> 
  987.       a myinvlaplace[ f, s, t, rm, rp, st, op ] /;
  988.     FreeQ[a, s],
  989.  
  990.   (*      -.  perform MyApart to redo partial fractions for    *)
  991.   (*          denominators with repeated roots that were not    *)
  992.   (*          properly handled by Apart.  kludge         *)
  993.   myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
  994.     partialFractionsKludge[x, s, t, rm, rp, st, op] /;
  995.     GetStateField[st, Lmyapartfield] && RationalPolynomialQ[x, s],
  996.  
  997.  
  998.   (*      D.  additivity                    *)
  999.   myinvlaplace[ x_ + y_, s_, t_, rm_, rp_, st_, op_ ] :>
  1000.       myinvlaplace[ x, s, t, rm, rp, st, op ] + 
  1001.     myinvlaplace[ y, s, t, rm, rp, st, op ], 
  1002.  
  1003.   myinvlaplace[ (x_ + y_) / c_, s_, t_, rm_, rp_, st_, op_ ] :>
  1004.       myinvlaplace[ x / c, s, t, rm, rp, st, op ] + 
  1005.     myinvlaplace[ y / c, s, t, rm, rp, st, op ], 
  1006.  
  1007.   myinvlaplace[ a_. Summation[i_, ib_, ie_, inc_][x_], s_, t_,
  1008.         rm_, rp_, st_, op_ ] :>
  1009.     Summation[i, ib, ie, inc][ myinvlaplace[ a x, s, t, rm, rp, st, op ] ],
  1010.  
  1011.  
  1012.   (*      E.  multiplication by an exponential            *)
  1013.   myinvlaplace[ a_. Exp[b_. s_ + c_.], s_, t_, rm_, rp_, st_, op_ ] :>
  1014.     substitutefort[ myinvlaplace[ a Exp[c], s, t, rm, rp, st, op ],
  1015.             t, t + b ] /;
  1016.     FreeQ[b, s],
  1017.  
  1018.  
  1019.   (*      F.  similarity                    *)
  1020.   (*          only use real, positive scaling factors        *)
  1021.   myinvlaplace[ f_, s_, t_, rm_, rp_, st_, op_ ] :>
  1022.     Block [ {newf, newrm, newrp, scale},
  1023.         scale = similarityScale[f, s];
  1024.         Assuming[ Positive[scale], op];
  1025.         newf = f /. s -> s / scale;
  1026.         newrm = If [ SameQ[rm, -Infinity], rm, scale rm ];
  1027.         newrp = If [ SameQ[rp, Infinity], rp, scale rp ];
  1028.         substitutefort[ myinvlaplace[newf, s, t, newrm, newrp, st, op],
  1029.                 t, t / scale ] / scale ] /;
  1030.     similarityQ[f, s],
  1031.  
  1032.  
  1033.   (*      G.  shift                        *)
  1034.   myinvlaplace[ f_, s_, t_, rm_, rp_, st_, op_ ] :>
  1035.     Block [    {newrm, newrp, normexpr, shift},
  1036.         {shift, normexpr} = GetShiftFactor[f, s];
  1037.         normexpr = normexpr /. s -> (s - shift);
  1038.         newrm = rm + Re[shift];
  1039.         newrp = rp + Re[shift];
  1040.         Exp[- shift t] *
  1041.           myinvlaplace[normexpr, s, t, newrm, newrp, st, op] ] /;
  1042.     ! SameQ[ First[GetShiftFactor[f, s]], 0 ],
  1043.  
  1044.  
  1045.   (*      H.  pole in denominator                *)
  1046.   myinvlaplace[ f_ / ( s_ + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
  1047.     myinvlaplace[ partialFractionsDialogue[f / (s + a), s, op], s, t,
  1048.               rm, rp, SetStateField[st, Lapartfield, False], op ] /;
  1049.     GetStateField[st, Lapartfield] && FreeQ[a, s],
  1050.  
  1051.   myinvlaplace[ f_ / ( s_ + a_ ), s_, t_, rm_, rp_, st_, op_ ] :>
  1052.     Exp[- a t] myinvlaplace[(f /. s -> s - a) / s, s, t, rm, rp, st, op] /;
  1053.     FreeQ[a, s],
  1054.  
  1055.  
  1056.   (*      I.  multiplication by s is progressive        *)
  1057.   (*          differentiation in the time domain;        *)
  1058.   (*          progressiveness handled by postinvrules        *)
  1059.   myinvlaplace[ s_^k_. a_., s_, t_, -Infinity, Infinity, st_, op_ ] :>
  1060.     Block [ {dialogue},
  1061.                 dialogue = SameQ[ Replace[Dialogue, op], All ];
  1062.         Assuming[ k > 0, dialogue ];
  1063.         If [ dialogue && ! IntegerQ[k],
  1064.              Print[ "where ", k, " is an integer" ] ];
  1065.           derivative[ myinvlaplace[a, s, t, rm, rp, st, op],
  1066.                 t, k, False ] ] /;
  1067.     FreeQ[k, s] && Implies[ NumberQ[k], IntegerQ[k] && ( k > 0 )],
  1068.  
  1069.   myinvlaplace[ s_^k_. a_., s_, t_, 0, Infinity, st_, op_ ] :>
  1070.     Block [ {dialogue},
  1071.                 dialogue = SameQ[ Replace[Dialogue, op], All ];
  1072.         Assuming[ k > 0, dialogue ];
  1073.         If [ dialogue && ! IntegerQ[k],
  1074.              Print[ "where ", k, " is an integer" ] ];
  1075.           derivative[myinvlaplace[a, s, t, rm, rp, st, op], t, k, 0] ] /;
  1076.     FreeQ[k, s] && Implies[ NumberQ[k], IntegerQ[k] && ( k > 0 )],
  1077.  
  1078.  
  1079.   (*      J.  division by s is integration in the time domain    *)
  1080.   myinvlaplace[ f_ s_^n_, s_, t_, rm_, rp_, st_, op_ ] :>
  1081.     Block [    {tau = Unique["tau"]},
  1082.         integrate[ myinvlaplace[f s^(n + 1), s, tau, rm, rp, st, op],
  1083.                t, tau ] ] /;
  1084.     IntegerQ[n] && (n < 0),
  1085.  
  1086.  
  1087.  
  1088.  
  1089.   (* V.   Strategies                        *)
  1090.  
  1091.   (*      A.  partial fractions                    *)
  1092.   myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
  1093.     myinvlaplace[ partialFractionsDialogue[x, s, op], s, t, rm, rp,
  1094.               SetStateField[st, Lapartfield, False], op ] /;
  1095.     GetStateField[st, Lapartfield],
  1096.  
  1097.  
  1098.   (*      B. Normalize denominator polynomials            *)
  1099.   myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
  1100.     myinvlaplace[ normalizeRatPoly[x, s], s, t, rm, rp,
  1101.               SetStateField[st, Lnormalizefield, False], op ] /;
  1102.     GetStateField[st, Lnormalizefield] && unnormalizedq[Denominator[x], s],
  1103.  
  1104.  
  1105.   (*      C.  factor denominator                *)
  1106.   myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
  1107.        Block [    {denom, numer},
  1108.         denom = Denominator[x];
  1109.         numer = Numerator[x];
  1110.         denom = Factor[denom];
  1111.         myinvlaplace[ numer / denom, s, t, rm, rp,
  1112.                   SetStateField[st, Lfactorfield, False], op ] ] /;
  1113.     GetStateField[st, Lfactorfield],
  1114.  
  1115.  
  1116.   (*      D.  distribute expression                *)
  1117.   myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
  1118.       myinvlaplace[ Distribute[x], s, t, rm, rp, st, op ] /;
  1119.     SameQ[Head[x], Times] && ! SameQ[x, Distribute[x]],
  1120.  
  1121.  
  1122.   (*      E.  expand numerators of expression            *)
  1123.   myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
  1124.       myinvlaplace[ Expand[x], s, t, rm, rp,
  1125.               SetStateField[st, Lexpandfield, False], op ] /;
  1126.     GetStateField[st, Lexpandfield],
  1127.  
  1128.  
  1129.   (*      F.  expand all terms in expression            *)
  1130.   myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
  1131.       myinvlaplace[ ExpandAll[x], s, t, rm, rp,
  1132.               SetStateField[st, Lexpandallfield, False], op ] /;
  1133.     GetStateField[st, Lexpandallfield],
  1134.  
  1135.  
  1136.   (*      G.  replace new operators and functions        *)
  1137.   (*          with their mathematical definition        *)
  1138.   myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
  1139.            myinvlaplace[ TheFunction[x], s, t, rm, rp,
  1140.               SetStateField[st, Lthefunfield, False], op ] /;
  1141.     GetStateField[st, Lthefunfield],
  1142.  
  1143.  
  1144.   (*      H.  apply definition if the user has enabled that    *)
  1145.   (*          the Definition option.                *)
  1146.   myinvlaplace[ x_, s_, t_, rm_, rp_, st_, op_ ] :>
  1147.     Block [    {state, trans},
  1148.         state = SetStateField[st, Ldefinitionfield, False];
  1149.         trans = Integrate[x Exp[s t], {s, rm, rp} ];
  1150.         If [ SameQ[Head[trans], Integrate],
  1151.              myinvlaplace[x, s, t, rm, rp, state, op],
  1152.              definitionDialogue[x, trans/(2 Pi), Integrate, op] ] ] /;
  1153.     GetStateField[st, Ldefinitionfield] && Replace[Definition, op]
  1154.  
  1155. }
  1156.  
  1157. (*  E N D     R U L E S  *)
  1158.  
  1159.  
  1160. (*  E N D     P A C K A G E  *)
  1161.  
  1162. End[]
  1163. EndPackage[]
  1164.  
  1165. If [ TrueQ[ $VersionNumber >= 2.0 ],
  1166.      On[ General::spell ];
  1167.      On[ General::spell1 ] ];
  1168.  
  1169.  
  1170. (*  A L I A S E S  *)
  1171.  
  1172. Unprotect[ { InverseLaPlace, InverseLaPlaceTransform } ]
  1173. InverseLaPlace = InvLaPlace
  1174. InverseLaPlaceTransform = InvLaPlace
  1175. Protect[ { InverseLaPlace, InverseLaPlaceTransform } ]
  1176.  
  1177.  
  1178. (*  H E L P     I N F O R M A T I O N  *)
  1179.  
  1180. Combine[ SPfunctions, { InvLaPlace } ]
  1181. Protect[ InvLaPlace ]
  1182.  
  1183.  
  1184. (*  E N D I N G     M E S S A G E  *)
  1185.  
  1186. Print[ "The inverse Laplace transform rule base InvLaPlace has been loaded." ]
  1187. Null
  1188.